CONTENTS
1. Introduction
Recent advances in next-generation sequencing, microarrays and mass spectrometry for omics data production have enabled the generation and collection of different modalities of high-dimensional molecular data1. Clustering multi-omic data has the potential to reveal further systems-level insights, but raises computational and biological challenges2. This vignette aims to show how to use the package named MOVICS to perform multi-omics integrative clustering and visualization for cancer subtyping researches. This R package provides a unified interface for 10 state-of-the-art multi-omics clustering algorithms, and standardizes the output for each algorithm so as to form a pipeline for downstream analyses. Ten algorithms are CIMLR, iClusterBayes, MoCluster, COCA, ConsensusClustering, IntNMF, LRAcluster, NEMO, PINSPlus, and SNF where the former three methods can also perform the process of feature selection. For cancer subtyping studies, MOVICS also integrates many downstream analyses for further subtype characterization and creates editable publication-quality illustrations (see more details below). Please note that MOVICS currently supports up to 6 omics data for jointly clustering and users must provide at least 2 omics datastes as input. Ok then, let us make our hands dirty to figure out how MOVICS works.
2. Installation
It is essential that you have R 3.6 or above already installed on your computer or server. MOVICS is a pipeline that utilises many other R packages that are currently available from CRAN and Bioconductor. For all of the steps of the pipeline to work, make sure that you have upgraded Bioconductor to newest version (BiocManager v3.10). After you have R and Bioconductor installed properly, you can start to install MOVICS. The easiest way to install it is by typing the following code into your R session:
When you are installing MOVICS, you may encounter some errors saying that some packages are not installed. These errors are caused by recursively depending on R packages, so if one package was not installed properly on your computer, MOVICS would fail. To solve these errors, you simply need to check those error messages, find out which packages required are missing, then install it with command BiocManager::install("YourErrorPackage") or install.packages("YourErrorPackage") directly. After that, retry installing MOVICS, it may take several times, but eventually it should work. After installation, you should be able to load the MOVICS package in your R session:
3. Real Data Scenario
This package contains two preprocessed real datasets of breast cancer. One dataset is brca.tcga.rda which is a list that includes 643 samples with four complete omics data types of breast cancer retrieved from TCGA-BRCA cohort3 (i.e., mRNA expression, lncRNA expression, DNA methylation profiling and somatic mutation matrix), and corresponding clinicopathological information (i.e., age, pathological stage, PAM50 subtype, vital status and overall suvival time); such data list also contains corresponding RNA-Seq raw count table and Fragments Per Kilobase Million (FPKM) data in order to test the functions for downstream analyses (e.g., differential expression analysis, drug sensitivity analysis, etc). The other one, brca.yau.rda, is an external validation dataset which contains gene expression profilings and clinicopathological information that downloaded from BRCA-YAU cohort4 with 682 samples (one sample without annotation of PAM50 subtype was removed), which can be used to test the predictive functions available in MOVICS. These two datasets can be loaded like below:
# load example data of breast cancer
load(system.file("extdata", "brca.tcga.rda", package = "MOVICS", mustWork = TRUE))
load(system.file("extdata", "brca.yau.rda", package = "MOVICS", mustWork = TRUE))Since in most cases, multi-omics clustering procedure is rather time-consuming, the TCGA-BRCA dataset was threfore first pre-processed to extract top 500 mRNAs, 500 lncRNA, 1,000 promoter CGI probes/genes with high variation using statistics of median absolute deviation (MAD), and 30 genes that mutated in at least 3% of the entire cohort. All these features are used in the following clustering analyses.
4. MOVICS Pipeline
4.1 Pipeline Introduction
MOVICS Pipeline diagram above outlines the concept for this package and such pipeline comprises separate functions that can be run individually (see datails about argument list by typing ?function_name() in R session). In above graphical abstract, all functions of MOVICS are included, which are categorized into three modules using three colors:
- GET Module: get subtypes through multi-omics integrative clustering
getElites(): get elites which are those features that pass the filtering procedure and are used for analysesgetClustNum(): get optimal cluster number by calculating clustering prediction index (CPI) and Gap-statisticsget%algorithm_name%(): get results from one specific multi-omics integrative clustering algorithm with detailed parametersgetMOIC(): get a list of results from multiple multi-omics integrative clustering algorithm with parameters by defaultgetConsensusMOIC(): get a consensus matrix that indicates the clustering robustness across different clustering algorithmsgetStdiz(): get a standardized data for generating comprehensive multi-omics heatmapgetMoHeatmap(): get a comprehensive multi-omics heatmap based on clustering results
- COMP Module: compare subtypes from multiple perspectives
compSurv(): compare survival outcome and generate a Kalan-Meier curve with pair-wise comparsion if possiblecompClinvar(): compare and summarize clinical features among different identified subtypescompMut(): compare mutational frequency and generate an oncoprint with significant mutationscompTMB(): compare total mutation burden among subtypes and generate distribution of Transitions and TransversionscompFGA(): compare fraction genome altered among subtypes and generate a barplot for distribution comparisoncompDrugsen(): compare estimated half maximal inhibitory concentration (\(IC_{50}\)) for drug sensitity and generate a boxviolin for distribution comparisoncompAgree(): compare agreement of current subtypes with other pre-existed classifications and generate an alluvial diagram and an agreement barplot
- RUN Module: run marker identification and verify subtypes
runDEA(): run differential expression analysis with three popular methods for choosing, including edgeR, DESeq2, and limmarunMarker(): run biomarker identification to determine uniquely and significanly differential expressed genes for each subtyperunGSEA(): run gene set enrichment analysis (GSEA), calculate activity of functional pathways and generate a pathway-specific heatmaprunNTP(): run nearest template prediction based on identified biomarkers to evaluate subtypes in external cohorts
4.2 Steps of Pipeline
Basically, the above three connected modules explain the workflow of this R package. MOVICS first identifies the cancer subtype (CS) by using one or multiple clustering algorithms; if multiple clustering algorithms are specified, it is highly recommended to perform a consensus clustering based on different subtyping results in order to derive stable and robust subtypes. Second, after having subtypes it is natural to exploit the heterogeneity of subtypes from as many angles as possible. Third, each subtype should have a list of subtype-specific biomarkers for reproducing such subtype in external cohorts. Therefore, let us follow this analytic pipeline and step into each module and each function in MOVICS.
4.2.1 GET Module
1) get data from example files
Data used for this vignette should be first extracted from the list of brca.tcga, including 4 types of multi-omics data. Except for omics data, this list also contains RNA-Seq raw count, FPKM matrix, MAF data, segmented copy number, clinical and survival information for downstream analyses. Notably, all these omics data used for clustering share exactly the same samples with exactly the same order, and please make sure of this when using MOVICS on your own data.
# print name of example data
names(brca.tcga)
#> [1] "mRNA.expr" "lncRNA.expr" "meth.beta" "mut.status" "count"
#> [6] "fpkm" "maf" "segment" "clin.info"
names(brca.yau)
#> [1] "mRNA.expr" "clin.info"
# extract multi-omics data
mo.data <- brca.tcga[1:4]
# extract raw count data for downstream analyses
count <- brca.tcga$count
# extract fpkm data for downstream analyses
fpkm <- brca.tcga$fpkm
# extract maf for downstream analysis
maf <- brca.tcga$maf
# extract segmented copy number for downstream analyses
segment <- brca.tcga$segment
# extract survival information
surv.info <- brca.tcga$clin.info2) get elites by reducing data dimension
Although all these omics data have been already processed (filtered from the original dataset), I am still pleased to show you how to use getElites() function to filter out features that meet some stringent requirements, and those features that are preserved in this procedure are considered elites by MOVICS. Four filtering methods are provided here, namely mad for median absolute deviation, sd for standard deviation, cox for univariate Cox proportional hazards regression, and none for nothing but is really necessary for binary omics data. This function also handles missing values coded in NA by removing them directly or imputing them by \(k\) nearest neighbors using a Euclidean metric through argument of na.action. Let me show you how to use getElites() below:
# scenario 1: considering we are dealing with an expression data that have 2 rows with NA values
tmp <- brca.tcga$mRNA.expr # get expression data
dim(tmp) # check data dimension
#> [1] 500 643
tmp[1,1] <- tmp[2,2] <- NA # set 2 rows with NA values
tmp[1:3,1:3] # check data
#> BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A
#> SCGB2A2 NA 1.42 7.24
#> SCGB1D2 10.11 NA 5.88
#> PIP 4.54 2.59 4.35
elite.tmp <- getElites(dat = tmp,
method = "mad",
na.action = "rm", # NA values will be removed
elite.pct = 1) # elite.pct equals to 1 means all (100%) features after NA removal will be selected even using mad method
#> --2 features with NA values are removed.
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # check dimension again and see that we have removed 2 rows with NA data
#> [1] 498 643
elite.tmp <- getElites(dat = tmp,
method = "mad",
na.action = "impute", # NA values will be imputed
elite.pct = 1)
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # all data kept
#> [1] 500 643
elite.tmp$elite.dat[1:3,1:3] # NA values have been imputed
#> BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A
#> SCGB2A2 6.867 1.420 7.24
#> SCGB1D2 10.110 4.739 5.88
#> PIP 4.540 2.590 4.35
# scenario 2: considering we are dealing with continous data and use mad or sd to select elites
tmp <- brca.tcga$mRNA.expr # get expression data with 500 features
elite.tmp <- getElites(dat = tmp,
method = "mad",
elite.pct = 0.1) # this time only top 10% features with high mad values are kept
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # get 50 elite left
#> [1] 50 643
elite.tmp <- getElites(dat = tmp,
method = "sd",
elite.num = 100, # this time only top 100 features with high sd values are kept
elite.pct = 0.1) # this time elite.pct arugument will be disabled because elite.num has been already indicated.
#> elite.num has been provided then discards elite.pct.
dim(elite.tmp$elite.dat) # get 100 elites left
#> [1] 100 643
# scenario 3: considering we are dealing with data and use cox to select elite
tmp <- brca.tcga$mRNA.expr # get expression data
elite.tmp <- getElites(dat = tmp,
method = "cox",
surv.info = surv.info, # must provide surivival information with 'futime' and 'fustat'
p.cutoff = 0.05,
elite.num = 100) # this time elite.num arugument will be disabled because cox method refers to p.cutoff to select elites
#> --all sample matched between omics matrix and survival data.
#> 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% 55% 60% 65% 70% 75% 80% 85% 90% 95% 100%
dim(elite.tmp$elite.dat) # get 125 elites
#> [1] 125 643
table(elite.tmp$unicox$pvalue < 0.05) # 125 genes have nominal pvalue < 0.05 in univariate Cox regression
#>
#> FALSE TRUE
#> 375 125
tmp <- brca.tcga$mut.status # get mutation data
elite.tmp <- getElites(dat = tmp,
method = "cox",
surv.info = surv.info, # must provide surivival information with 'futime' and 'fustat'
p.cutoff = 0.05,
elite.num = 100) # this time elite.num arugument will be disabled because cox method refers to p.cutoff to select elites
#> --all sample matched between omics matrix and survival data.
#> 7% 13% 20% 27% 33% 40% 47% 53% 60% 67% 73% 80% 87% 93% 100%
dim(elite.tmp$elite.dat) # get 3 elites
#> [1] 3 643
table(elite.tmp$unicox$pvalue < 0.05) # 3 mutations have nominal pvalue < 0.05
#>
#> FALSE TRUE
#> 27 3
# scenario 4: considering we are dealing with mutation data without using cox to select elites
tmp <- brca.tcga$mut.status # get mutation data
rowSums(tmp) # check mutation frequency
#> PIK3CA TP53 TTN CDH1 GATA3 MLL3 MUC16 MAP3K1 SYNE1 MUC12
#> 208 186 111 83 58 49 48 38 33 32
#> DMD NCOR1 FLG PTEN RYR2 USH2A SPTA1 MAP2K4 MUC5B NEB
#> 31 31 30 29 27 27 25 25 24 24
#> SPEN MACF1 RYR3 DST HUWE1 HMCN1 CSMD1 OBSCN APOB SYNE2
#> 23 23 23 22 22 22 21 21 21 21
elite.tmp <- getElites(dat = tmp,
method = "none", # must set as 'none'
elite.num = 80, # note: in this scenario elite.num refer to frequency of mutation
elite.pct = 0.1) # discard because elite.num has been already indicated
#> --method of none only supports binary omics data (e.g., somatic mutation matrix), and in this manner, elite.pct and elite.num are used to cut frequency.
#> elite.num has been provided then discards elite.pct.
rowSums(elite.tmp$elite.dat) # only genes that are mutated in over than 80 samples are kept as elites
#> PIK3CA TP53 TTN CDH1
#> 208 186 111 83
elite.tmp <- getElites(dat = tmp,
method = "none", # must set as 'none'
elite.pct = 0.2) # note: in this scenario elite.pct refer to frequency of mutation / sample size
#> --method of none only supports binary omics data (e.g., somatic mutation matrix), and in this manner, elite.pct and elite.num are used to cut frequency.
#> missing elite.num then use elite.pct
rowSums(elite.tmp$elite.dat) # only genes that are mutated in over than 0.2*643=128.6 samples are kept as elites
#> PIK3CA TP53
#> 208 186
# get mo.data list just like below (not run)
# mo.data <- list(omics1 = elite.tmp$elite.dat,
# omics2 = ...)Now I think I have made this clear for you concerning how to reduce data dimension for clustering analysis. Since the mo.data has been already prepared, I am going to stick on this data list, and take you further to MOVICS.
3) get optimal number for clustering
The most important parameter to estimate in any clustering study is the optimum number of clusters \(k\) for the data, where \(k\) needs to be small enough to reduce noise but large enough to retain important information. Herein MOVICS refers to CPI5 and Gaps-statistics6 to estimate the number of clusters by using getClustNum() function.
# identify optimal clustering number (may take a while)
optk.brca <- getClustNum(data = mo.data,
is.binary = c(F,F,F,T), # note: the 4th data is somatic mutation which is a binary matrix
try.N.clust = 2:8, # try cluster number from 2 to 8
fig.name = "CLUSTER NUMBER OF TCGA-BRCA")
#> calculating Cluster Prediction Index...
#> 5% complete
#> 5% complete
#> 10% complete
#> 10% complete
#> 15% complete
#> 15% complete
#> 20% complete
#> 25% complete
#> 25% complete
#> 30% complete
#> 30% complete
#> 35% complete
#> 35% complete
#> 40% complete
#> 45% complete
#> 45% complete
#> 50% complete
#> 50% complete
#> 55% complete
#> 55% complete
#> 60% complete
#> 65% complete
#> 65% complete
#> 70% complete
#> 70% complete
#> 75% complete
#> 75% complete
#> 80% complete
#> 85% complete
#> 85% complete
#> 90% complete
#> 90% complete
#> 95% complete
#> 95% complete
#> 100% complete
#> calculating Gap-statistics...
#> visualization done...
#> --the imputed optimal cluster number is 3 arbitrarily, but it would be better referring to other priori knowledge.
Figure 1. Identification of optimal cluster number by calculating Cluster Prediction Index (blue line) and Gaps-statistics (red line) in TCGA-BRCA cohort.
The above estimation of clustering number gives an arbitrary \(k\) of 2. However, the popular PAM50 classifier for breast cancer has 5 classifications and if taking a close look at the descriptive figure, it can be noticed that both CPI and Gaps-statistics does not decline too much at \(k\) of 5. Therefore, considering these knowledge, \(k\) of 5 is chosen as the optimal clustering number for futher analyses.
4) get results from single algorithm
In this part I will first show how to use MOVICS to perform multi-omics integrative clustering by specifying one algorithm with detailed parameters. For example, let us try iClusterBayes method like below:
# perform iClusterBayes (may take a while)
iClusterBayes.res <- getiClusterBayes(data = mo.data,
N.clust = 5,
type = c("gaussian","gaussian","gaussian","binomial"),
n.burnin = 1800,
n.draw = 1200,
prior.gamma = c(0.5, 0.5, 0.5, 0.5),
sdev = 0.05,
thin = 3)
#> clustering done...
#> feature selection done...Otherwise, a unified function can be used for all algorithms individually with detailed parameters, that is, getMOIC().
iClusterBayes.res <- getMOIC(data = mo.data,
N.clust = 5,
methodslist = "iClusterBayes", # specify only ONE algorithm here
type = c("gaussian","gaussian","gaussian","binomial"), # data type corresponding to the list of matrics
n.burnin = 1800,
n.draw = 1200,
prior.gamma = c(0.5, 0.5, 0.5, 0.5),
sdev = 0.05,
thin = 3)By specifying only one algorithm (i.e., iClusterBayes) in the argument of methodslist, the above getMOIC() will return exactly the same results from getiClusterBayes() if the same parameters are provided. The returned result contains a clust.res object that has two columns: clust to indicate the subtype that the sample belongs to, and samID records the corresponding sample name. For algorithms that provide feature selection procedure (i.e., iClusterBayes, CIMLR, and MoCluster), the result also contains a feat.res object that stores the information of such procedure. To those algorithms involving hierarchical clustering (e.g., COCA, ConsensusClustering), the corresponding dendrogram for sample clustering will be also returned as clust.dend.
5) get results from multiple algorithms at once
If you simultaneously specify a list of algorithms to methodslist argument in getMOIC(), it will automatically perform each algorithm with default parameters one by one, and a list of results derived from specified algorithms will be finally returned. Now that iClusterBayes has been finished, let us try other 9 algorithms at once with parameters by default. This will take some time, so have a coffee break.
# perform multi-omics integrative clustering with the rest of 9 algorithms
moic.res.list <- getMOIC(data = mo.data,
methodslist = list("SNF", "PINSPlus", "NEMO", "COCA", "LRAcluster", "ConsensusClustering", "IntNMF", "CIMLR", "MoCluster"),
N.clust = 5,
type = c("gaussian", "gaussian", "gaussian", "binomial"))
#>
|
| | 0%
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%
# attach iClusterBayes.res as a list using append() to moic.res.list with 9 results already
moic.res.list <- append(moic.res.list,
list("iClusterBayes" = iClusterBayes.res))
# save moic.res.list
save(moic.res.list, file = "moic.res.list.rda")6) get consensus from different algorithms
Since now all the clustering results from 10 algorthms are in hand, MOVICS borrows the idea of consensus clustering for later integration of the clustering results derived from different algorithms, so as to improve the clustering robustness. To be specific, if \(t_{max}\) algorithms are specified where \(2\le t_{max} \le10\), getConsensusMOIC() calculates a matrix \(M_{n\times n}^{(t)}\) per algorithm where \(n\) is the number of samples, and \(M_{ij}^{(t)}=1\) only when the sample \(i\) and \(j\) are clustered in the same subtype, otherwise \(M_{ij}^{(t)}=0\). After get all results from specified algorithms, MOVICS calculates a consensus matrix \(CM=\sum_{t=1}^{t_{max}}M^{(t)}\), and \(cm_{ij}\in[0,10]\). Applying hierarchical clustering to \(CM\), MOVICS generates a clustering result which represents a robust pair-wise similarities for samples because it considers different multi-omics integrative clustering algorithms. The simplest way to perform getConsensusMOIC() is to pass a list of object returned by multiple get%algorithm_name%() or by getMOIC() with specific argument of methodslist. This can be done like below:
cmoic.brca <- getConsensusMOIC(moic.res.list = moic.res.list,
fig.name = "CONSENSUS HEATMAP",
distance = "euclidean",
linkage = "average")
Figure 2. Consensus heatmap based on results from 10 multi-omics integrative clustering algorithms with cluster number of 5.
Remember, you must choose at least two methods to perform the above consensus clustering by getConsensusMOIC(), and this function will return a consensus matrix (a probability matrix represents how many times samples belonging to the same subtype can be clustered together by different multi-omics clustering methods) and a corresponding consensus heatmap. Ideally, the consensus heatmap will show a perfect diagonal rectangle, and the input values are 0 and 1 only because all algorithms derived the same clustering results.
7) get multi-omics heatmap based on clustering result
Genome-wide heatmaps are widely used to graphically display potential underlying patterns within the large genomic dataset. They have been used to reveal information about how the samples/genes cluster together and provide insights into potential sample biases or other artifacts. Herein MOVICS provides getMoHeatmap to visually deal with preclustered multi-omics data and generate an exquisite heatmap for publication requirement. Before using getMoHeatmap(), omics data should be properly processed by using the function of getStdiz() which returns a list storing normalized omics data. Omics data, especially for expression (e.g., RNA and protein), should be centered (centerFlag = T) or scaled (scaleFlag = T) or z-scored (both centered and scaled). Generally, DNA methylation data (\(\beta\) matrix ranging from 0 to 1) and somatic mutation (0 and 1 binary matrix) should not be normalized. However, it is a good choice to convert the methylation \(\beta\) value to M value following the formula of \(M=log_2\frac{\beta}{1-\beta}\) for its stronger signal in visualization and M value is suitable for normalization. This function also provides an argument of halfwidth for continuous omics data; such argument is used to truncate the ‘extremum’ after normalization; specifically, normalized values that exceed the halfwidth boundaries will be replaced by the halfwidth, which is vary useful to map colors in heatmap.
# convert beta value to M value for stronger signal
indata <- mo.data
indata$meth.beta <- log2(indata$meth.beta/(1 - indata$meth.beta))
# data normalization for heatmap
plotdata <- getStdiz(data = indata,
halfwidth = c(2,2,2,NA), # no truncation for mutation
centerFlag = c(T,T,T,F), # no center for mutation
scaleFlag = c(T,T,T,F)) # no scale for mutationAs I mentioned earlier, several algorithms also provide feature selection; those selected features show a complex cross-talk with other omics data and might have special biological significance that drives the heterogeneity of cancers. Therefore I show below how to generate a comprehensive heatmap based on a single algorithm (e.g., iClusterBayes) with selected features by using getMoHeatmap(). However in the first place, features must be extracted:
feat <- iClusterBayes.res$iClusterBayes$feat.res
feat1 <- feat[which(feat$dataset == "mRNA.expr"),][1:10,"feature"]
feat2 <- feat[which(feat$dataset == "lncRNA.expr"),][1:10,"feature"]
feat3 <- feat[which(feat$dataset == "meth.beta"),][1:10,"feature"]
feat4 <- feat[which(feat$dataset == "mut.status"),][1:10,"feature"]
annRow <- list(feat1, feat2, feat3, feat4)The feat.res contained in iClusterBayes.res is sorted by posterior probability of features for each omics data. In this manner, the top 10 features are selected for each omics data and a feature list is generated and named as annRow for heatmap raw annotation (Do not worry, I will talk about how to attach column annotation for samples in a minute).
# set color for each omics data
# if no color list specified all subheatmaps will be unified to green and red color pattern
mRNA.col <- c("#00FF00", "#008000", "#000000", "#800000", "#FF0000")
lncRNA.col <- c("#6699CC", "white" , "#FF3C38")
meth.col <- c("#0074FE", "#96EBF9", "#FEE900", "#F00003")
mut.col <- c("grey90" , "black")
col.list <- list(mRNA.col, lncRNA.col, meth.col, mut.col)
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = c("mRNA","lncRNA","Methylation","Mutation"),
is.binary = c(F,F,F,T), # the 4th data is mutation which is binary
legend.name = c("mRNA.FPKM","lncRNA.FPKM","M value","Mutated"),
clust.res = iClusterBayes.res$clust.res, # cluster results
clust.dend = NULL, # no dendrogram
show.rownames = c(F,F,F,F), # specify for each omics data
show.colnames = F, # show no sample names
annRow = annRow, # mark selected features
color = col.list,
annCol = NULL, # no annotation for samples
annColors = NULL, # no annotation color
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF ICLUSTERBAYES")
Figure 3. Comprehensive heatmap of multi-omics integrative clustering by iClusterBayes algorithm with annotation of potential cluster drivers (selected features).
Of course, since there are a total of 10 results stored in the moic.res.list, you can also choose any one of them to create the heatmap. Here I select COCA which also returns sample dendrogram like below:
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = c("mRNA","lncRNA","Methylation","Mutation"),
is.binary = c(F,F,F,T), # the 4th data is mutation which is binary
legend.name = c("mRNA.FPKM","lncRNA.FPKM","M value","Mutated"),
clust.res = moic.res.list$COCA$clust.res, # cluster results
clust.dend = moic.res.list$COCA$clust.dend, # show dendrogram for samples
color = col.list,
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF COCA")
Figure 4. Comprehensive heatmap of multi-omics integrative clustering by COCA algorithm with dendrogram for samples.
Now go back to the consensus result of cmoic.brca that integrates 10 algorithms and this time samples’ annotation is also provided to generate the heatmap. Since the core function of getMoHeatmap() is based on ComplexHeatmap R package, when creating annotations, you should always use circlize::colorRamp2() function to generate the color mapping function for continuous variables (e.g., age in this example).
# extract PAM50, pathologic stage and age for sample annotation
annCol <- surv.info[,c("PAM50", "pstage", "age"), drop = F]
# generate corresponding colors for sample annotation
annColors <- list(age = circlize::colorRamp2(breaks = c(min(annCol$age),
median(annCol$age),
max(annCol$age)),
colors = c("#0000AA", "#555555", "#AAAA00")),
PAM50 = c("Basal" = "blue",
"Her2" = "red",
"LumA" = "yellow",
"LumB" = "green",
"Normal" = "black"),
pstage = c("T1" = "green",
"T2" = "blue",
"T3" = "red",
"T4" = "yellow",
"TX" = "black"))
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = c("mRNA","lncRNA","Methylation","Mutation"),
is.binary = c(F,F,F,T), # the 4th data is mutation which is binary
legend.name = c("mRNA.FPKM","lncRNA.FPKM","M value","Mutated"),
clust.res = cmoic.brca$clust.res, # consensusMOIC results
clust.dend = NULL, # show no dendrogram for samples
show.rownames = c(F,F,F,F), # specify for each omics data
show.colnames = F, # show no sample names
show.row.dend = c(F,F,F,F), # show no dendrogram for features
annRow = NULL, # no selected features
color = col.list,
annCol = annCol, # annotation for samples
annColors = annColors, # annotation color
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF CONSENSUSMOIC")
Figure 5. Comprehensive heatmap of multi-omics integrative clustering based on consensus across 10 algorithms with clinicopathological annotation for samples.
4.2.2 COMP Module
After identification of cancer subtypes, it is essential to further characterize each subtype by discovering their difference from multile aspects. To this end, MOVICS provides commonly used downstream analyses in cancer subtyping researches for easily cohesion with results derived from GET Module. Now let us check them out.
1) compare survival outcome
I start with comparing the prognosis among different subtypes. MOVICS provides function of compSurv() which not only calculates the overall nominal P value by log-rank test, but also performs pair-wise comparsion and derives adjusted P values if more than two subtypes are identified. These information will be all printed in the Kaplan-Meier Curve which is convenient for researchers to refer. Except for clustering results (e.g., cmoic.brca in this example), you must additionally provide surv.info argument which should be a data.frame that stores a futime column for survival time (unit of day) and another fustat column for survival outcome (0 = censor; 1 = event).
# survival comparsion
surv.brca <- compSurv(moic.res = cmoic.brca,
surv.info = surv.info,
convt.time = "m", # convert day unit to month
surv.median.line = "h",
fig.name = "KAPLAN-MEIER CURVE OF CONSENSUSMOIC") # draw horizontal line at median survival
#> --a total of 643 samples are identified.
#> --removed missing values.
#> --leaving 642 observations.
#> --cut survival curve up to 10 years.
Figure 6. Kaplan-Meier survival curve of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
print(surv.brca)
#> $fitd
#> Call:
#> survdiff(formula = Surv(futime, fustat) ~ Subtype, data = mosurv.res,
#> na.action = na.exclude)
#>
#> N Observed Expected (O-E)^2/E (O-E)^2/V
#> Subtype=CS1 146 28 14.6 12.269 15.455
#> Subtype=CS2 132 17 15.1 0.251 0.317
#> Subtype=CS3 107 7 13.4 3.028 3.756
#> Subtype=CS4 144 7 17.4 6.218 8.175
#> Subtype=CS5 113 16 14.6 0.140 0.176
#>
#> Chisq= 22.3 on 4 degrees of freedom, p= 2e-04
#>
#> $fit
#> Call: survfit(formula = Surv(futime, fustat) ~ Subtype, data = mosurv.res,
#> na.action = na.exclude, error = "greenwood", type = "kaplan-meier",
#> conf.type = "plain")
#>
#> n events median 0.95LCL 0.95UCL
#> CS1 146 28 130 67.3 NA
#> CS2 132 17 114 83.6 144
#> CS3 107 7 NA 102.5 NA
#> CS4 144 7 216 113.5 NA
#> CS5 113 16 NA 97.2 NA
#>
#> $overall.p
#> [1] 0.000172211
#>
#> $pairwise.p
#>
#> Pairwise comparisons using Log-Rank test
#>
#> data: mosurv.res and Subtype
#>
#> CS1 CS2 CS3 CS4
#> CS2 0.11909 - - -
#> CS3 0.01345 0.11909 - -
#> CS4 0.00055 0.03497 0.72167 -
#> CS5 0.16073 0.87433 0.16073 0.03022
#>
#> P value adjustment method: BH2) compare clinical features
Then compare the clinical features among different subtypes. MOVICS provides function of compClinvar() which enables summarizing both continous and categorical variables and performing proper statistical tests. This function can give a table in .docx format that is easy to use in medical research papers.
clin.brca <- compClinvar(moic.res = cmoic.brca,
var2comp = surv.info, # data.frame needs to summarize
strata = "Subtype", # stratifying (subtype) variable
factorVars = c("PAM50","pstage","fustat"), # features that are considered categorical variables
nonnormalVars = "futime", # feature(s) that are considered using nonparametric test
exactVars = "pstage", # feature(s) that are considered using exact test
doWord = T, # generate .docx file
tab.name = "SUMMARIZATION OF CLINICAL FEATURES")
#> --all samples matched.
#> Warning in jstable::CreateTableOne2(vars = setdiff(colnames(dat),
#> strata), : NAs introduced by coercion| level | CS1 | CS2 | CS3 | CS4 | CS5 | p | test | sig | |
|---|---|---|---|---|---|---|---|---|---|
| n | 146 | 133 | 107 | 144 | 113 | ||||
| fustat (%) | 0 | 118 (80.8) | 115 (86.5) | 100 (93.5) | 137 (95.1) | 97 (85.8) | 0.001 | ** | |
| 1 | 28 (19.2) | 18 (13.5) | 7 ( 6.5) | 7 ( 4.9) | 16 (14.2) | ||||
| futime (median [IQR]) | 719.00 [431.50, 1355.00] | 746.50 [432.25, 1420.25] | 703.00 [436.00, 1270.50] | 848.50 [470.50, 1547.75] | 742.00 [470.00, 1605.00] | 0.691 | nonnorm | ||
| PAM50 (%) | Basal | 2 ( 1.4) | 0 ( 0.0) | 0 ( 0.0) | 0 ( 0.0) | 109 (96.5) | <0.001 | ** | |
| Her2 | 38 (26.0) | 0 ( 0.0) | 0 ( 0.0) | 0 ( 0.0) | 0 ( 0.0) | ||||
| LumA | 47 (32.2) | 78 (58.6) | 91 (85.0) | 128 (88.9) | 0 ( 0.0) | ||||
| LumB | 53 (36.3) | 55 (41.4) | 16 (15.0) | 0 ( 0.0) | 0 ( 0.0) | ||||
| Normal | 6 ( 4.1) | 0 ( 0.0) | 0 ( 0.0) | 16 (11.1) | 4 ( 3.5) | ||||
| pstage (%) | T1 | 40 (27.4) | 35 (26.3) | 32 (29.9) | 40 (27.8) | 25 (22.1) | NA | exact | |
| T2 | 85 (58.2) | 77 (57.9) | 63 (58.9) | 72 (50.0) | 70 (61.9) | ||||
| T3 | 12 ( 8.2) | 17 (12.8) | 10 ( 9.3) | 31 (21.5) | 14 (12.4) | ||||
| T4 | 9 ( 6.2) | 4 ( 3.0) | 2 ( 1.9) | 1 ( 0.7) | 3 ( 2.7) | ||||
| TX | 0 ( 0.0) | 0 ( 0.0) | 0 ( 0.0) | 0 ( 0.0) | 1 ( 0.9) | ||||
| age | 58.47 (13.75) | 59.32 (14.03) | 59.62 (12.40) | 57.58 (12.23) | 55.67 (12.11) | 0.116 |
3) compare mutational frequency
Subype-specific mutation might be promising as therapeutic target. Hence, I then compare the mutational frequency among different clusters. MOVICS provides function of compMut() which deals with binary mutational data (0 = wild; 1 = mutated). This function applies independent testing (i.e., Fisher' s exact test or \(\chi^2\) test) for each mutation, and generates a table with statistical results (.docx file also if specified), and creates an oncoprint using those differentially mutated genes.
# mutational frequency comparison
mut.brca <- compMut(moic.res = cmoic.brca,
mut.matrix = brca.tcga$mut.status, # binary somatic mutation matrix
doWord = T, # generate table in .docx format
doPlot = T, # draw oncoprint
freq.cutoff = 0.05, # keep those genes that mutated in at least 5% of samples
p.adj.cutoff = 0.05, # keep those genes with adjusted p value < 0.05 to draw oncoprint
innerclust = T, # perform clustering within each subtype
annCol = annCol, # same annotation for heatmap
annColors = annColors, # same annotation color for heatmap
width = 6,
height = 2,
fig.name = "ONCOPRINT FOR SIGNIFICANT MUTATIONS",
tab.name = "INDEPENDENT TEST BETWEEN SUBTYPE AND MUTATION")
#> --all samples matched.
Figure 7. Mutational oncoprint of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
| Gene (Mutated) | TMB | CS1 | CS2 | CS3 | CS4 | CS5 | pvalue | padj |
|---|---|---|---|---|---|---|---|---|
| PIK3CA | 208 (32%) | 55 (37.7%) | 10 ( 7.5%) | 77 (72.0%) | 64 (44.4%) | 2 ( 1.8%) | 6.90e-42 | 3.10e-41 |
| TP53 | 186 (29%) | 83 (56.8%) | 12 ( 9.0%) | 0 ( 0.0%) | 8 ( 5.6%) | 83 (73.5%) | 2.06e-63 | 1.85e-62 |
| TTN | 111 (17%) | 41 (28.1%) | 14 (10.5%) | 14 (13.1%) | 22 (15.3%) | 20 (17.7%) | 2.03e-03 | 3.04e-03 |
| CDH1 | 83 (13%) | 11 ( 7.5%) | 3 ( 2.3%) | 19 (17.8%) | 49 (34.0%) | 1 ( 0.9%) | 5.38e-19 | 1.61e-18 |
| GATA3 | 58 ( 9%) | 5 ( 3.4%) | 32 (24.1%) | 10 ( 9.3%) | 11 ( 7.6%) | 0 ( 0.0%) | 6.21e-11 | 1.40e-10 |
| MLL3 | 49 ( 8%) | 11 (7.5%) | 12 (9.0%) | 10 (9.3%) | 11 (7.6%) | 5 (4.4%) | 6.31e-01 | 6.31e-01 |
| MUC16 | 48 ( 8%) | 18 (12.3%) | 8 ( 6.0%) | 6 ( 5.6%) | 8 ( 5.6%) | 8 ( 7.1%) | 2.05e-01 | 2.31e-01 |
| MAP3K1 | 38 ( 6%) | 2 ( 1.4%) | 6 ( 4.5%) | 16 (15.0%) | 12 ( 8.3%) | 2 ( 1.8%) | 3.65e-05 | 6.57e-05 |
| SYNE1 | 33 ( 5%) | 9 (6.2%) | 1 (0.8%) | 5 (4.7%) | 8 (5.6%) | 10 (8.8%) | 3.27e-02 | 4.20e-02 |
4) compare total mutation burden
Needless to say, immunotherapy is becoming a pillar of modern cancer treatment. Recent analyses have linked the tumoral genomic landscape with antitumor immunity. In particular, it has been proposed that overall mutational load might drive T-cell responses7, whereas tumor aneuploidy correlates with markers of immune evasion and reduced response to immunotherapy8. To quantify these genomic alterations that may affect immunotherapy, MOVICS provides two functions to calculate total mutation burden (TMB) and fraction genome altered (FGA). To be specific, Total Mutations Burden (TMB) refers to the number of mutations that are found in the tumor genome, while FGA is the percentage of genome that has been affected by copy number gains or losses . Both attributes are useful for genetic researchers as they provide them with more in-depth information on the genomic make-up of the tumors. Let me show you how to use these two functions by starting with compTMB(). First of all, the input maf data for this function must have the following 10 columns at least:
head(maf)
#> Tumor_Sample_Barcode Hugo_Symbol Chromosome Start_Position End_Position
#> 1 BRCA-A1XY-01A USP24 chr1 55159655 55159655
#> 2 BRCA-A1XY-01A ERICH3 chr1 74571494 74571494
#> 3 BRCA-A1XY-01A KIF26B chr1 245419680 245419680
#> 4 BRCA-A1XY-01A USP34 chr2 61189055 61189055
#> 5 BRCA-A1XY-01A ANTXR1 chr2 69245305 69245305
#> 6 BRCA-A1XY-01A SCN9A chr2 166199365 166199365
#> Variant_Classification Variant_Type Reference_Allele Tumor_Seq_Allele1
#> 1 Missense_Mutation SNP T T
#> 2 Missense_Mutation SNP C C
#> 3 Silent SNP G G
#> 4 Silent SNP G G
#> 5 Silent SNP G G
#> 6 Silent SNP G G
#> Tumor_Seq_Allele2
#> 1 C
#> 2 T
#> 3 T
#> 4 C
#> 5 A
#> 6 ARun compTMB() then. By default, this function considers nonsynonymous variants only when counting somatic mutation frequency, including Frame_Shift_Del, Frame_Shift_Ins, Splice_Site, Translation_Start_Site, Nonsense_Mutation, Nonstop_Mutation, In_Frame_Del, In_Frame_Ins, and Missense_Mutation by default. In addition to TMB, this function also classifies Single Nucleotide Variants into Transitions and Transversions (TiTv), and depicts the distribution of both TMB and TiTv.
# compare TMB
tmb.brca <- compTMB(moic.res = cmoic.brca,
maf = maf,
rmDup = T, # remove duplicated variants per sample
rmFLAGS = F, # keep FLAGS mutations
exome.size = 38, # estimated exome size
test.method = "nonparametric", # statistical testing method
fig.name = "DISTRIBUTION OF TMB AND TITV")
#> --67 samples mismatched from current subtypes.
#> -Validating
#> -Silent variants: 24329
#> -Summarizing
#> --Possible FLAGS among top ten genes:
#> TTN
#> MUC16
#> -Processing clinical data
#> --Missing clinical data
#> -Finished in 4.350s elapsed (3.710s cpu)
Figure 8. Comparison of TMB and TiTv among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
5) compare fraction genome altered
Next, compFGA() calculates not only FGA but also computes specific gain (FGG) or loss (FGL) per sample within each subtype. To get this function worked, an eligible input of segmented copy number should be prepared with exactly the same column name like below:
# change column names of segment data
colnames(segment) <- c("sample","chrom","start","end","value")Let’s see how the input looks like:
head(segment)
#> sample chrom start end value
#> 1 BRCA-A090-01A 1 3301765 54730235 -0.1271
#> 2 BRCA-A090-01A 1 54730247 57443819 -0.0899
#> 3 BRCA-A090-01A 1 57448465 57448876 -1.1956
#> 4 BRCA-A090-01A 1 57448951 64426102 -0.1009
#> 5 BRCA-A090-01A 1 64426648 106657734 -0.1252
#> 6 BRCA-A090-01A 1 106657854 106835667 0.1371Then this function can be run without any difficulties. Notably, if your CNA calling procedure did not provide a segmented copy number as value column but the original copy number, argument of iscopynumber must be switched to TRUE instead.
# compare FGA, FGG, and FGL
fga.brca <- compFGA(moic.res = cmoic.brca,
segment = segment,
iscopynumber = F, # this is a segmented copy number file
cnathreshold = 0.2, # threshold to determine CNA gain or loss
test.method = "nonparametric", # statistical testing method
fig.name = "BARPLOT OF FGA")
#> --2 samples mismatched from current subtypes.
#> 5% 10% 15% 21% 26% 31% 36% 41% 46% 51% 57% 62% 67% 72% 77% 82% 88% 93% 98%
Figure 9. Barplot of fraction genome altered among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
| samID | FGA | FGG | FGL | Subtype |
|---|---|---|---|---|
| BRCA-A03L-01A | 0.6217991 | 0.3086727 | 0.3131264 | CS1 |
| BRCA-A04R-01A | 0.2531019 | 0.0913201 | 0.1617818 | CS2 |
| BRCA-A075-01A | 0.7007067 | 0.4144424 | 0.2862643 | CS1 |
| BRCA-A08O-01A | 0.6501287 | 0.4564814 | 0.1936473 | CS3 |
| BRCA-A0A6-01A | 0.1468893 | 0.0635649 | 0.0833244 | CS4 |
| BRCA-A0AD-01A | 0.1722214 | 0.0386452 | 0.1335762 | CS2 |
6) compare drug sensitivity
Predicting response to medication is particularly important for drugs with a narrow therapeutic index, for example chemotherapeutic agents, because response is highly variable and side effects are potentially lethal. Therefore, Paul Geeleher et al. (2014)9 used baseline gene expression and \(in\;vitro\) drug sensitivity derived from cell lines, coupled with \(in\;vivo\) baseline tumor gene expression, to predict patients’ response to drugs. Paul developped an R package pRRophetic for prediction of clinical chemotherapeutic response from tumor gene expression levels10, and now this function has been involved in MOVICS to examine difference of drug sensitivity among different clusters. Here I estimate the \(IC_{50}\) of Cisplatin for 5 identified subtypes of breast cancer in TCGA.
# drug sensitivity comparsion
drug.brca <- compDrugsen(moic.res = cmoic.brca,
norm.expr = fpkm[,cmoic.brca$clust.res$samID], # double guarantee sample order
drugs = c("Cisplatin", "Paclitaxel"), # a vector of names of drug in GDSC
tissueType = "breast", # choose specific tissue type to construct ridge regression model
test.method = "nonparametric", # statistical testing method
prefix = "BOXVIOLIN OF ESTIMATED IC50")
#> --all samples matched.
#> --log2 transformation done for expression data.
#> Cisplatin done...
#> Paclitaxel done...
Figure 10. Boxviolins for estimated IC50 of Cisplatin and Paclitaxel among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
| Est.IC50 | Subtype | |
|---|---|---|
| BRCA-A03L-01A | 3.801326 | CS1 |
| BRCA-A075-01A | 3.834292 | CS1 |
| BRCA-A0AW-01A | 3.127552 | CS1 |
| BRCA-A0AZ-01A | 3.302353 | CS1 |
| BRCA-A0BC-01A | 3.409796 | CS1 |
| BRCA-A0BF-01A | 3.357585 | CS1 |
7) compare agreement with other subtypes
At present, many cancers have traditional classifications, and evaluating the consistency of new subtypes with previous classifications is critical to reflect the robustness of clustering analysis and to determine potential but novel subtypes. To measure the agreement (similarity) between the current subtypes and other pre-existed classifications MOVICS provides function of compAgree() to calculate four statistics: Rand Index (RI)11, Adjusted Mutual Information (AMI)12, Jaccard Index (JI)13, and Fowlkes-Mallows (FM)14; all these measurements range from 0 to 1 and the larger the value is, the more similar the two appraises are. This function can also generate an alluvial diagram to visualize the agreement of two appraises with the current subtypes as reference.
# customize the factor level for pstage
surv.info$pstage <- factor(surv.info$pstage, levels = c("TX","T1","T2","T3","T4"))
# agreement comparsion (support up to 6 classifications include current subtype)
agree.brca <- compAgree(moic.res = cmoic.brca,
subt2comp = surv.info[,c("PAM50","pstage")],
doPlot = T,
box.width = 0.2,
fig.name = "AGREEMENT OF CONSENSUSMOIC WITH PAM50 AND PSTAGE")
#> --all samples matched.
Figure 11. Agreement of 5 identified subtypes of breast cancer with PAM50 classification and pathological stage in TCGA-BRCA cohort.
| current.subtype | other.subtype | RI | AMI | JI | FM |
|---|---|---|---|---|---|
| Subtype | PAM50 | 0.6929744 | 0.4043456 | 0.2910887 | 0.4694389 |
| Subtype | pstage | 0.5506751 | 0.0124902 | 0.1565996 | 0.2884961 |
4.2.3 RUN Module
Take a deep breath, we successfully enter the last RUN Module, and victory is in sight. In this module, MOVICS aims to characterize different subtypes by identifying their potential predictive biomarkers and functional pathways. Identifying and applying molecular biomarkers to predict subtype with efficiency is particularly important for disease management and treatment, thus improving clinical outcome. In this context, MOVICS searches for subtype-specific biomarkers by starting with differential expression analysis (DEA).
1) run differential expression analysis
MOVICS provides runDEA() function which embeds three state-of-the-art DEA approches to identify differentially expressed genes (DEGs), including edgeR and DESeq2 for RNA-Seq count data and limma for microarray data. Notably, since runDEA() checks the data scale automatically when choosing limma algorithm, it is recommended to provide a normalized expression data (e.g., RSEM, FPKM, TPM) witout z-score or log2 transformation. This step is also quite time-consuming, especially in the case of using raw count data, so release your hand from the keyboard, close your eyes, and relax for a moment.
# run DEA with edgeR
runDEA(dea.method = "edger",
expr = count, # raw count data
moic.res = cmoic.brca,
prefix = "TCGA-BRCA") # prefix of figure name
#> --all samples matched.
#> --you choose edger and please make sure an RNA-Seq count data was provided.
#> edger of CS1_vs_Others done...
#> edger of CS2_vs_Others done...
#> edger of CS3_vs_Others done...
#> edger of CS4_vs_Others done...
#> edger of CS5_vs_Others done...
# run DEA with DESeq2
runDEA(dea.method = "deseq2",
expr = count,
moic.res = cmoic.brca,
prefix = "TCGA-BRCA")
#> --all samples matched.
#> --you choose deseq2 and please make sure an RNA-Seq count data was provided.
#> deseq2 of CS1_vs_Others done...
#> deseq2 of CS2_vs_Others done...
#> deseq2 of CS3_vs_Others done...
#> deseq2 of CS4_vs_Others done...
#> deseq2 of CS5_vs_Others done...
# run DEA with limma
runDEA(dea.method = "limma",
expr = fpkm, # normalized expression data
moic.res = cmoic.brca,
prefix = "TCGA-BRCA")
#> --all samples matched.
#> --you choose limma and please make sure a normalized expression data [FPKM or TPM without log2 transformation is recommended] was provided.
#> --log2 transformation done for expression data.
#> limma of CS1_vs_Others done...
#> limma of CS2_vs_Others done...
#> limma of CS3_vs_Others done...
#> limma of CS4_vs_Others done...
#> limma of CS5_vs_Others done...Each identified cancer subtype will be compared with the rest (Others) and the corresponding .txt file will be stored according to argument of res.path. By default, these files will be saved under the current working directory.
2) run biomarker identification procedure
In this procedure, the most differentially expressed genes sorted by log2FoldChange are chosen as the biomarkers for each subtype. These biomarkers should pass the significance threshold (e.g., nominal \(P\) value < 0.05 and adjusted \(P\) value < 0.05) and must not overlap with any biomarkers indentified for other subtypes.
# choose edgeR result to identify subtype-specific up-regulated biomarkers
marker.up <- runMarker(moic.res = cmoic.brca,
dea.method = "edger", # name of DEA method
prefix = "TCGA-BRCA", # MUST be the same of argument in runDEA()
dat.path = getwd(), # path of DEA files
res.path = getwd(), # path to save marker files
p.cutoff = 0.05, # p cutoff to identify significant DEGs
p.adj.cutoff = 0.05, # padj cutoff to identify significant DEGs
dirct = "up", # direction of dysregulation in expression
n.marker = 100, # number of biomarkers for each subtype
doplot = T, # generate diagonal heatmap
norm.expr = fpkm, # use normalized expression as heatmap input
annCol = annCol, # sample annotation in heatmap
annColors = annColors, # colors for sample annotation
show_rownames = F, # show no rownames (biomarker name)
fig.name = "UPREGULATED BIOMARKER HEATMAP")
#> --all samples matched.
#> --log2 transformation done for expression data.
Figure 12. Heatmap of subtype-specific upregulated biomarkers using edgeR algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
| probe | class | dirct |
|---|---|---|
| PNMT | CS1 | up |
| AKR1B15 | CS1 | up |
| DLK1 | CS1 | up |
| ACE2 | CS1 | up |
| CRISP3 | CS1 | up |
| ACSM1 | CS1 | up |
Then try results derived from limma to identify subtype-specific downregulated biomarkers. Please feel free to try DESeq2.
# choose DESeq2 result to identify subtype-specific down-regulated biomarkers
marker.dn <- runMarker(moic.res = cmoic.brca,
dea.method = "limma",
prefix = "TCGA-BRCA",
dirct = "down",
n.marker = 50, # switch to 50
doplot = T,
norm.expr = fpkm,
annCol = annCol,
annColors = annColors,
fig.name = "DOWNREGULATED BIOMARKER HEATMAP")
#> --all samples matched.
#> --log2 transformation done for expression data.
Figure 13. Heatmap of subtype-specific downregulated biomarkers using limma algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
3) run gene set enrichment analysis
Similarly, GSEA is run for each subtype based on its corresponding DEA result to identify subtype-specific functional pathways. To this end, I have prepared a geneset background which includes all gene sets derived from GO biological processes (c5.bp.v7.1.symbols.gmt) from The Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).
# MUST locate ABSOLUTE path of msigdb file
MSIGDB.FILE <- system.file("extdata", "c5.bp.v7.1.symbols.xls", package = "MOVICS", mustWork = TRUE)Likewise, these identified specific pathways should pass the significance threshold (e.g., nominal \(P\) value < 0.05 and adjusted \(P\) value < 0.25) and must not overlap with any pathways identified for other subtypes. After having the subtype-specific pathways, genes that are inside the pathways are retrieved to calculate a single sample enrichment score by using GSVA R package. Subsequently, subtype-specific enrichment score will be represented by the mean or median (mean by default) value within the subtype, and will be further visualized by diagonal heatmap.
# run GSEA to identify up-regulated GO pathways using results from edgeR
gsea.up <- runGSEA(moic.res = cmoic.brca,
dea.method = "edger", # name of DEA method
prefix = "TCGA-BRCA", # MUST be the same of argument in runDEA()
dat.path = getwd(), # path of DEA files
res.path = getwd(), # path to save GSEA files
msigdb.path = MSIGDB.FILE, # MUST be the ABSOLUTE path of msigdb file
norm.expr = fpkm, # use normalized expression to calculate enrichment score
dirct = "up", # direction of dysregulation in pathway
p.cutoff = 0.05, # p cutoff to identify significant pathways
p.adj.cutoff = 0.25, # padj cutoff to identify significant pathways
gsva.method = "gsva", # method to calculate single sample enrichment score
norm.method = "mean", # normalization method to calculate subtype-specific enrichment score
fig.name = "UPREGULATED PATHWAY HEATMAP")
#> --all samples matched.
#> Loading required package: GSEABase
#> Loading required package: annotate
#> Loading required package: AnnotationDbi
#> Loading required package: stats4
#> Loading required package: IRanges
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:dplyr':
#>
#> first, rename
#> The following object is masked from 'package:base':
#>
#> expand.grid
#>
#> Attaching package: 'IRanges'
#> The following objects are masked from 'package:dplyr':
#>
#> collapse, desc, slice
#> The following object is masked from 'package:grDevices':
#>
#> windows
#>
#> Attaching package: 'AnnotationDbi'
#> The following object is masked from 'package:dplyr':
#>
#> select
#> Loading required package: XML
#> Loading required package: graph
#>
#> Attaching package: 'graph'
#> The following object is masked from 'package:XML':
#>
#> addNode
#> GSEA done...
#> --log2 transformation done for expression data.
#> Estimating GSVA scores for 50 gene sets.
#> Computing observed enrichment scores
#> Estimating ECDFs with Gaussian kernels
#> Using parallel with 1 cores
#>
|
| | 0%
|
|= | 2%
|
|=== | 4%
|
|==== | 6%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 14%
|
|========== | 16%
|
|============ | 18%
|
|============= | 20%
|
|============== | 22%
|
|================ | 24%
|
|================= | 26%
|
|================== | 28%
|
|==================== | 30%
|
|===================== | 32%
|
|====================== | 34%
|
|======================= | 36%
|
|========================= | 38%
|
|========================== | 40%
|
|=========================== | 42%
|
|============================= | 44%
|
|============================== | 46%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 54%
|
|==================================== | 56%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================== | 64%
|
|=========================================== | 66%
|
|============================================ | 68%
|
|============================================== | 70%
|
|=============================================== | 72%
|
|================================================ | 74%
|
|================================================= | 76%
|
|=================================================== | 78%
|
|==================================================== | 80%
|
|===================================================== | 82%
|
|======================================================= | 84%
|
|======================================================== | 86%
|
|========================================================= | 88%
|
|========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 94%
|
|============================================================== | 96%
|
|================================================================ | 98%
|
|=================================================================| 100%
#> gsva done...
#> heatmap done...
Figure 14. Heatmap of subtype-specific upregulated pathways using edgeR algorithm for 5 identified subtypes in TCGA-BRCA cohort.
Check some columns of GSEA results for the first subtype (CS1).
| setSize | enrichmentScore | NES | pvalue | |
|---|---|---|---|---|
| GO_REGULATION_OF_ION_TRANSPORT | 420 | -0.4832370 | -1.555051 | 0.0010040 |
| GO_EMBRYONIC_MORPHOGENESIS | 459 | -0.4504464 | -1.449967 | 0.0010050 |
| GO_EPITHELIAL_CELL_DIFFERENTIATION | 470 | -0.4956895 | -1.597430 | 0.0010050 |
| GO_SYNAPTIC_SIGNALING | 441 | -0.5216787 | -1.677297 | 0.0010050 |
| GO_BEHAVIOR | 392 | -0.5060951 | -1.625033 | 0.0010091 |
| GO_REGULATION_OF_TRANS_SYNAPTIC_SIGNALING | 293 | -0.5122532 | -1.631490 | 0.0010101 |
Also check results of subtype-specific enrichment scores.
| CS1 | CS2 | CS3 | CS4 | CS5 | |
|---|---|---|---|---|---|
| GO_BENZENE_CONTAINING_COMPOUND_METABOLIC_PROCESS | 0 | 0 | 0 | 0 | -1 |
| GO_REGULATION_OF_MONOCYTE_CHEMOTAXIS | 0 | 0 | 0 | 0 | 0 |
Then try results derived from DESeq2 to identify subtype-specific downregulated pathways. Please feel free to try limma.
# run GSEA to identify down-regulated GO pathways using results from DESeq2
gsea.dn <- runGSEA(moic.res = cmoic.brca,
dea.method = "deseq2",
prefix = "TCGA-BRCA",
msigdb.path = MSIGDB.FILE,
norm.expr = fpkm,
dirct = "down",
p.cutoff = 0.05,
p.adj.cutoff = 0.25,
gsva.method = "ssgsea", # switch to ssgsea
norm.method = "median", # switch to median
fig.name = "DOWNREGULATED PATHWAY HEATMAP")
#> --all samples matched.
#> GSEA done...
#> --log2 transformation done for expression data.
#> Estimating ssGSEA scores for 50 gene sets.
#>
|
| | 0%Using parallel with 1 cores
#>
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======= | 12%
|
|======== | 12%
|
|======== | 13%
|
|========= | 13%
|
|========= | 14%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 19%
|
|============= | 20%
|
|============= | 21%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|=================== | 30%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|====================== | 35%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 36%
|
|======================== | 37%
|
|======================== | 38%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 45%
|
|============================== | 46%
|
|============================== | 47%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 53%
|
|=================================== | 54%
|
|=================================== | 55%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================= | 61%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 62%
|
|========================================= | 63%
|
|========================================= | 64%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 70%
|
|============================================== | 71%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|================================================== | 78%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 79%
|
|==================================================== | 80%
|
|==================================================== | 81%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|======================================================== | 87%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================= | 95%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 96%
|
|=============================================================== | 97%
|
|=============================================================== | 98%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%
#> ssgsea done...
#> heatmap done...
Figure 15. Heatmap of subtype-specific downregulated pathways using limma algorithm for 5 identified subtypes in TCGA-BRCA cohort.
4) run nearest template prediction in external cohort
Oh wait, did we forget something? Yeah there is a dataset that has not been used yet, so let’s see if we can use those subtype-specific biomarkers to verify the current breast cancer subtypes in the external Yau cohort. In this part, our core purpose is to predict the possible subtypes of each sample in the external dataset. In most cases, this is a multi-classification problem, and the identified biomarkers may be difficult to be overall matched in the external cohort, so it might not be reliable to use model-based predictive algorithms. Therefore, MOVICS switches to nearest template prediction (NTP) which can be flexibly applied to cross-platform, cross-species, and multiclass predictions without any optimization of analysis parameters15. Only one thing have to do is generating a template file, which has been fortunately prepared already.
# run NTP in Yau cohort by using up-regulated biomarkers
brca.pred <- runNTP(expr = brca.yau$mRNA.expr,
templates = marker.up$templates, # the template has been already prepared in runMarker()
scale = T, # scale input data
center = T, # center input data
doPlot = T, # to generate heatmap
fig.name = "NTP HEATMAP FOR YAU")
#> --original template has 500 biomarkers and 272 are matched in external expression profiling.
#> cosine correlation distance
#> 682 samples; 5 classes; 39-66 features/class
#> serial processing; 1000 permutation(s)...
#> predicted samples/class (FDR<0.05)
Figure 16. Heatmap of NTP in Yau cohort using subtype-specific upregulated biomarkers identified from TCGA-BRCA cohort
#>
#> CS1 CS2 CS3 CS4 CS5 <NA>
#> 104 85 90 120 140 143| prediction | d.CS1 | d.CS2 | d.CS3 | d.CS4 | d.CS5 | p.value | FDR | |
|---|---|---|---|---|---|---|---|---|
| 107 | CS2 | 0.7344 | 0.5165 | 0.7377 | 0.7471 | 0.7512 | 0.001 | 0.0020 |
| 109 | CS1 | 0.5816 | 0.7581 | 0.7172 | 0.7308 | 0.7295 | 0.001 | 0.0020 |
| 11 | CS1 | 0.6527 | 0.7307 | 0.7430 | 0.7557 | 0.7731 | 0.001 | 0.0020 |
| 110 | CS2 | 0.7619 | 0.5418 | 0.7583 | 0.7655 | 0.7567 | 0.001 | 0.0020 |
| 111 | CS1 | 0.6806 | 0.7157 | 0.7198 | 0.7889 | 0.7911 | 0.007 | 0.0106 |
| 113 | CS4 | 0.6785 | 0.7505 | 0.7394 | 0.6389 | 0.7308 | 0.007 | 0.0106 |
Above an object brca.pred is prepared which is similar in structure of object returned by getMOIC() but only stores clust.res which can be passed to functions within COMP Module if additional data is available. For example, I herein first compare the survival outcome of the predicted 5 cancer subtypes in Yau cohort.
# compare survival outcome in Yau cohort
surv.yau <- compSurv(moic.res = brca.pred,
surv.info = brca.yau$clin.info,
convt.time = "m", # switch to year
surv.median.line = "hv", # switch to both
fig.name = "KAPLAN-MEIER CURVE OF NTP FOR YAU")
#> --a total of 682 samples are identified.
#> --cut survival curve up to 10 years.
Figure 17. Kaplan-Meier survival curve of predicted 5 subtypes of breast cancer in Yau cohort.
print(surv.yau)
#> $fitd
#> Call:
#> survdiff(formula = Surv(futime, fustat) ~ Subtype, data = mosurv.res,
#> na.action = na.exclude)
#>
#> N Observed Expected (O-E)^2/E (O-E)^2/V
#> Subtype=CS1 136 51 43.8 1.169 1.452
#> Subtype=CS2 117 51 38.5 4.037 4.876
#> Subtype=CS3 119 33 42.3 2.045 2.517
#> Subtype=CS4 159 43 57.3 3.590 4.820
#> Subtype=CS5 151 50 46.0 0.351 0.441
#>
#> Chisq= 11.2 on 4 degrees of freedom, p= 0.02
#>
#> $fit
#> Call: survfit(formula = Surv(futime, fustat) ~ Subtype, data = mosurv.res,
#> na.action = na.exclude, error = "greenwood", type = "kaplan-meier",
#> conf.type = "plain")
#>
#> n events median 0.95LCL 0.95UCL
#> CS1 136 51 NA NA NA
#> CS2 117 51 205 109 NA
#> CS3 119 33 236 NA NA
#> CS4 159 43 222 191 NA
#> CS5 151 50 NA NA NA
#>
#> $overall.p
#> [1] 0.02410534
#>
#> $pairwise.p
#>
#> Pairwise comparisons using Log-Rank test
#>
#> data: mosurv.res and Subtype
#>
#> CS1 CS2 CS3 CS4
#> CS2 0.699 - - -
#> CS3 0.158 0.059 - -
#> CS4 0.100 0.039 0.901 -
#> CS5 0.804 0.504 0.244 0.158
#>
#> P value adjustment method: BHI further check the agreement between the predicted subtype and PAM50 classification.
# compare agreemeng in Yau cohort
agree.yau <- compAgree(moic.res = brca.pred,
subt2comp = brca.yau$clin.info[, "PAM50", drop = F],
doPlot = T,
fig.name = "YAU PREDICTEDMOIC WITH PAM50")
#> --all samples matched.
Figure 18. Agreement of predicted 5 subtypes of breast cancer with PAM50 classification in Yau cohort.
| current.subtype | other.subtype | RI | AMI | JI | FM |
|---|---|---|---|---|---|
| Subtype | PAM50 | 0.7830213 | 0.3864876 | 0.3318526 | 0.4994425 |
It is clearly that the predicted subtypes of breast cancer in Yau cohort can distinguish prognosis well and also show similar pattern of agreement with PAM50 classification as compared to TCGA-BRCA to some extent.
5. Little Trick
Indeed, the moic.res parameter is necessary for almost all the functions, especially for downstream analyses. You may deem that if the moic.res object is not obtained properly through the module of GET, downstream analyses cannot run smoothly. But in fact, you can fool downstream modules (COMP and RUN) with little tricks. The core information stored in moic.res is a data.frame named clust.res with a samID column (character) for samples name (should be exactly the same with rownames) and a clust column (integer) for subtype indicator. Another information requried for some functions in downstream analyses is mo.method, a stringe value that is usually considered as prefix for output files. For example, if I want to check the prognostic value of PAM50 classification, the only thing I have to prepare is a pseudo moic.res object with a mo.method just like below:
# include original clinical information as `clust.res` and a string value for `mo.method` to a list
pseudo.moic.res <- list("clust.res" = surv.info,
"mo.method" = "PAM50")
# make pseudo samID
pseudo.moic.res$clust.res$samID <- rownames(pseudo.moic.res$clust.res)
# make pseudo clust using a mapping relationship
pseudo.moic.res$clust.res$clust <- sapply(pseudo.moic.res$clust.res$PAM50,
switch,
"Basal" = 1, # relabel Basal as 1
"Her2" = 2, # relabel Her2 as 2
"LumA" = 3, # relabel LumA as 3
"LumB" = 4, # relabel LumnB as 4
"Normal" = 5) # relabel Normal as 5Let check how the pseudo moic.res object looks like.
head(pseudo.moic.res$clust.res)
#> fustat futime PAM50 pstage age samID clust
#> BRCA-A03L-01A 0 2442 LumA T3 34 BRCA-A03L-01A 3
#> BRCA-A04R-01A 0 2499 LumB T1 36 BRCA-A04R-01A 4
#> BRCA-A075-01A 0 518 LumB T2 42 BRCA-A075-01A 4
#> BRCA-A08O-01A 0 943 LumA T2 45 BRCA-A08O-01A 3
#> BRCA-A0A6-01A 0 640 LumA T2 64 BRCA-A0A6-01A 3
#> BRCA-A0AD-01A 0 1157 LumA T1 83 BRCA-A0AD-01A 3Ok everything is ready, just keep in mind how your interested subtypes are mapped and feel free to use such pseudo object to perform downstream analyses such as compSurv() as follows:
# survival comparsion
pam50.brca <- compSurv(moic.res = pseudo.moic.res,
surv.info = surv.info,
convt.time = "y", # convert day unit to year
surv.median.line = "h", # draw horizontal line at median survival
fig.name = "KAPLAN-MEIER CURVE OF PAM50 BY PSEUDO")
#> --a total of 643 samples are identified.
#> --removed missing values.
#> --leaving 642 observations.
#> --cut survival curve up to 10 years.
Figure 19. Kaplan-Meier survival curve of PAM50 subtypes of breast cancer with pseudo input in TCGA-BRCA cohort.
6. Summary
I have introduced the pipeline and described the functions included in MOVICS in almost detail. This package is rather young and I hope to improve it further in the future, along with more functionality. If you have any questions, bug reports, or suggestions for improving MOVICS, please email them to xlu.cpu@foxmail.com. Feeback is crucial for improving a package that tries to seamlessly incorporate functionality and flexibility from many other useful tools.
7. Session Information
sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18362)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=Chinese (Simplified)_China.936
#> [3] LC_MONETARY=Chinese (Simplified)_China.936
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=Chinese (Simplified)_China.936
#>
#> attached base packages:
#> [1] stats4 parallel stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] GSEABase_1.48.0 graph_1.64.0 annotate_1.64.0
#> [4] XML_3.98-1.20 AnnotationDbi_1.48.0 IRanges_2.20.0
#> [7] S4Vectors_0.24.0 MOVICS_1.0.0 testthat_2.2.1
#> [10] Biobase_2.46.0 BiocGenerics_0.32.0 devtools_2.2.1
#> [13] usethis_1.5.1 kableExtra_1.1.2 knitr_1.25
#> [16] dplyr_0.8.3
#>
#> loaded via a namespace (and not attached):
#> [1] Hmisc_4.2-0 corpcor_1.6.9
#> [3] aricode_1.0.0 class_7.3-15
#> [5] ps_1.3.0 foreach_1.4.7
#> [7] rprojroot_1.3-2 crayon_1.3.4
#> [9] MASS_7.3-51.6 nlme_3.1-140
#> [11] backports_1.1.5 sva_3.34.0
#> [13] impute_1.60.0 GOSemSim_2.12.0
#> [15] rlang_0.4.7 readxl_1.3.1
#> [17] svd_0.5 XVector_0.26.0
#> [19] heatmap.plus_1.3 irlba_2.3.3
#> [21] nloptr_1.2.1 callr_3.3.2
#> [23] limma_3.42.0 BiocParallel_1.20.0
#> [25] rjson_0.2.20 bit64_0.9-7
#> [27] glue_1.3.1 rngtools_1.4
#> [29] SNFtool_2.3.0 processx_3.4.1
#> [31] NEMO_0.1.0 oompaData_3.1.1
#> [33] DOSE_3.12.0 haven_2.3.1
#> [35] tidyselect_0.2.5 SummarizedExperiment_1.16.0
#> [37] km.ci_0.5-2 rio_0.5.16
#> [39] tidyr_1.0.0 zoo_1.8-8
#> [41] ggpubr_0.2.3 maftools_2.2.0
#> [43] ridge_2.4 xtable_1.8-4
#> [45] magrittr_1.5 prettyGraphs_2.1.6
#> [47] evaluate_0.14 bibtex_0.4.2
#> [49] ggplot2_3.3.2 cli_1.1.0
#> [51] zlibbioc_1.32.0 rstudioapi_0.10
#> [53] rpart_4.1-15 wordcloud_2.6
#> [55] fastmatch_1.1-0 ClassDiscovery_3.3.12
#> [57] InterSIM_2.2.0 shiny_1.4.0
#> [59] GSVA_1.34.0 prettydoc_0.4.0
#> [61] xfun_0.10 coca_1.1.0
#> [63] clue_0.3-57 pkgbuild_1.0.6
#> [65] cluster_2.1.0 caTools_1.17.1.2
#> [67] tidygraph_1.1.2 tibble_3.0.3
#> [69] flexclust_1.4-0 ggrepel_0.8.1
#> [71] png_0.1-7 permute_0.9-5
#> [73] withr_2.1.2 bitops_1.0-6
#> [75] ggforce_0.3.1 cellranger_1.1.0
#> [77] plyr_1.8.6 e1071_1.7-2
#> [79] survey_4.0 pillar_1.4.6
#> [81] RcppParallel_4.4.4 geepack_1.3-1
#> [83] gplots_3.0.1.1 GlobalOptions_0.1.1
#> [85] fs_1.3.1 iClusterPlus_1.22.0
#> [87] GetoptLong_0.1.7 graphite_1.32.0
#> [89] europepmc_0.3 clusterProfiler_3.14.0
#> [91] vctrs_0.3.2 ellipsis_0.3.0
#> [93] generics_0.0.2 urltools_1.7.3
#> [95] NMF_0.21.0 tools_3.6.1
#> [97] foreign_0.8-71 entropy_1.2.1
#> [99] munsell_0.5.0 tweenr_1.0.1
#> [101] fgsea_1.12.0 DelayedArray_0.12.0
#> [103] abind_1.4-5 fastmap_1.0.1
#> [105] compiler_3.6.1 pkgload_1.0.2
#> [107] httpuv_1.5.2 sessioninfo_1.1.1
#> [109] pkgmaker_0.27 GenomeInfoDbData_1.2.2
#> [111] gridExtra_2.3 edgeR_3.28.0
#> [113] lattice_0.20-38 jstable_0.9.5
#> [115] later_1.0.0 alluvial_0.1-2
#> [117] jsonlite_1.6 scales_1.1.1
#> [119] dendsort_0.3.3 carData_3.0-2
#> [121] pbapply_1.4-2 genefilter_1.68.0
#> [123] promises_1.1.0 car_3.0-4
#> [125] SimDesign_2.0.1 doParallel_1.0.15
#> [127] latticeExtra_0.6-28 checkmate_1.9.4
#> [129] openxlsx_4.1.2 rmarkdown_2.3
#> [131] cowplot_1.0.0 statmod_1.4.32
#> [133] webshot_0.5.1 forcats_0.4.0
#> [135] igraph_1.2.4.1 survival_3.2-3
#> [137] IntNMF_1.2.0 yaml_2.2.0
#> [139] htmltools_0.4.0 memoise_1.1.0
#> [141] modeltools_0.2-22 locfit_1.5-9.1
#> [143] graphlayouts_0.5.0 viridisLite_0.3.0
#> [145] digest_0.6.22 assertthat_0.2.1
#> [147] mime_0.7 rappdirs_0.3.1
#> [149] registry_0.5-1 LRAcluster_1.0
#> [151] KMsurv_0.1-5 RSQLite_2.1.2
#> [153] remotes_2.1.0 data.table_1.12.6
#> [155] blob_1.2.0 vegan_2.5-6
#> [157] mogsa_1.20.0 CMScaller_0.99.1
#> [159] preprocessCore_1.48.0 survMisc_0.5.5
#> [161] labeling_0.3 shinythemes_1.1.2
#> [163] splines_3.6.1 Formula_1.2-3
#> [165] RCurl_1.95-4.12 broom_0.7.0
#> [167] hms_0.5.2 colorspace_1.4-1
#> [169] ConsensusClusterPlus_1.50.0 base64enc_0.1-3
#> [171] BiocManager_1.30.9 GenomicRanges_1.38.0
#> [173] shape_1.4.4 nnet_7.3-12
#> [175] Rcpp_1.0.2 mclust_5.4.5
#> [177] circlize_0.4.8 enrichplot_1.6.0
#> [179] R6_2.4.0 survminer_0.4.6
#> [181] grid_3.6.1 ggridges_0.5.1
#> [183] lifecycle_0.2.0 acepack_1.4.1
#> [185] labelled_2.5.0 oompaBase_3.2.9
#> [187] zip_2.0.4 curl_4.2
#> [189] ggsignif_0.6.0 minqa_1.2.4
#> [191] gdata_2.18.0 PINSPlus_2.0.4
#> [193] DO.db_2.9 Matrix_1.2-17
#> [195] qvalue_2.18.0 desc_1.2.0
#> [197] RColorBrewer_1.1-2 iterators_1.0.12
#> [199] stringr_1.4.0 htmlwidgets_1.5.1
#> [201] officer_0.3.12 polyclip_1.10-0
#> [203] triebeard_0.3.0 purrr_0.3.3
#> [205] CIMLR_1.0.0 gridGraphics_0.4-1
#> [207] rvest_0.3.4 ComplexHeatmap_2.5.4
#> [209] mgcv_1.8-28 patchwork_0.0.1
#> [211] htmlTable_1.13.2 bdsmatrix_1.3-4
#> [213] codetools_0.2-16 matrixStats_0.55.0
#> [215] GO.db_3.10.0 FNN_1.1.3
#> [217] gtools_3.8.1 prettyunits_1.0.2
#> [219] gridBase_0.4-7 GenomeInfoDb_1.22.0
#> [221] gtable_0.3.0 DBI_1.0.0
#> [223] ggpmisc_0.3.5 ggalluvial_0.12.1
#> [225] highr_0.8 httr_1.4.1
#> [227] KernSmooth_2.23-15 stringi_1.4.3
#> [229] progress_1.2.2 reshape2_1.4.3
#> [231] farver_2.0.3 uuid_0.1-4
#> [233] viridis_0.5.1 magick_2.2
#> [235] coxme_2.2-16 xml2_1.2.2
#> [237] rvcheck_0.1.6 boot_1.3-22
#> [239] lme4_1.1-21 geneplotter_1.64.0
#> [241] ggplotify_0.0.4 DESeq2_1.26.0
#> [243] bit_1.1-14 ExPosition_2.8.23
#> [245] ggraph_2.0.0 pkgconfig_2.0.3
#> [247] mitools_2.4 tableone_0.12.08. Citing MOVICS
MOVICS is a wrapper incorporating algorithms and functions from many other sources. The majority of functions used in MOVICS are incorporated directly from other packages or are drawn from specific published algorithms. Hence, below I provide guidance on which papers should be cited alongside MOVICS when using these functions:
If you use MOVICS R package in published research, please cite:
- HAVE NOT BEEN PUBLISHED
Please cite corresponding literature below for multi-omic clustering algorithm if used:
- CIMLR: Ramazzotti D, Lal A, Wang B, Batzoglou S, Sidow A (2018). Multi-omic tumor data reveal diversity of molecular mechanisms that correlate with survival. Nat Commun, 9(1):4453.
- iClusterBayes: Mo Q, Shen R, Guo C, Vannucci M, Chan KS, Hilsenbeck SG (2018). A fully Bayesian latent variable model for integrative clustering analysis of multi-type omics data. Biostatistics, 19(1):71-86.
- Mocluster: Meng C, Helm D, Frejno M, Kuster B (2016). moCluster: Identifying Joint Patterns Across Multiple Omics Data Sets. J Proteome Res, 15(3):755-765.
- COCA: Hoadley KA, Yau C, Wolf DM, et al (2014). Multiplatform analysis of 12 cancer types reveals molecular classification within and across tissues of origin. Cell, 158(4):929-944.
- ConsensusClustering: Monti S, Tamayo P, Mesirov J, et al (2003). Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data. Mach Learn, 52:91-118.
- IntNMF: Chalise P, Fridley BL (2017). Integrative clustering of multi-level omic data based on non-negative matrix factorization algorithm. PLoS One, 12(5):e0176278.
- LRAcluster: Wu D, Wang D, Zhang MQ, Gu J (2015). Fast dimension reduction and integrative clustering of multi-omics data using low-rank approximation: application to cancer molecular classification. BMC Genomics, 16(1):1022.
- NEMO: Rappoport N, Shamir R (2019). NEMO: cancer subtyping by integration of partial multi-omic data. Bioinformatics, 35(18):3348-3356.
- PINSPlus: Nguyen H, Shrestha S, Draghici S, Nguyen T (2019). PINSPlus: a tool for tumor subtype discovery in integrated genomic data. Bioinformatics, 35(16):2843-2846.
- SNF`: Wang B, Mezlini AM, Demir F, et al (2014). Similarity network fusion for aggregating data types on a genomic scale. Nat Methods, 11(3):333-337.
Heatmap generated by MOVICS uses ComplexHeatmap R package, so please cite if any heatmap is created:
- Gu Z, Eils R, Schlesner M (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 32(18):2847–2849.
If FLAGS is removed when calculating TMB, please cite the following two at the same time, otherwise only the first one:
- Mayakonda A, Lin D C, Assenov Y, et al. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res, 28(11):1747-1756.
- Shyr C, Tarailo-Graovac M, Gottlieb M, Lee JJ, van Karnebeek C, Wasserman WW. (2014). FLAGS, frequently mutated genes in public exomes. BMC Med Genomics, 7(1): 1-14.
If calculating FGA, please cite:
- Cerami E, Gao J, Dogrusoz U, et al. (2012). The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov, 2(5):401-404.
- Gao J, Aksoy B A, Dogrusoz U, et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal, 6(269):pl1-pl1.
Drug sensitivity analysis is based on pRRophetic R package, so please cite:
- Geeleher P, Cox N, Huang R S (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One, 9(9):e107468.
- Geeleher P, Cox N J, Huang R S (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol, 15(3):1-12.
Please cite corresponding literature below if using differential expression analysis:
- edgeR:
- Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1):139-140.
- McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40(10):4288-4297.
- DESeq2: Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15(12):550-558.
- limma: Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 43(7):e47.
Gene set enrichment analysis uses clusterProfiler R package with the following paper:
- Yu G, Wang L, Han Y, He Q (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5):284-287.
For single sample enrichment analysis, plese cite literatures below appropriately:
- ssgsea: Barbie, D.A. et al (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature, 462(5):108-112.
- gsva: Hänzelmann, S., Castelo, R. and Guinney, J. (2013). GSVA: Gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics, 14(1):7.
- zscore: Lee, E. et al (2008). Inferring pathway activity toward precise disease classification. PLoS Comp Biol, 4(11):e1000217.
- plage: Tomfohr, J. et al (2005). Pathway level analysis of gene expression using singular value decomposition. BMC Bioinformatics, 6(1):1-11.
External validation uses nearest template prediction as implanted in CMScaller R package. Please cite:
- Eide P W, Bruun J, Lothe R A, et al. (2017). CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models. Sci Rep, 7(1):1-8.
- Hoshida, Y. (2010). Nearest Template Prediction: A Single-Sample-Based Flexible Class Prediction with Confidence Assessment. PLoS One, 5(11):e15543.
REFERENCES
- Pierre-Jean M, Deleuze J F, Le Floch E, et al (2019). Clustering and variable selection evaluation of 13 unsupervised methods for multi-omics data integration. Brief Bioinformatics.
- Rappoport N, Shamir R (2018). Multi-omic and multi-view clustering algorithms: review and cancer benchmark. Nucleic Acids Res, 46(20):10546-10562.
- Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature, 2012, 490(7418):61-78.
- Yau C, Esserman L, Moore D H, et al. (2010). A multigene predictor of metastatic outcome in early stage hormone receptor-negative and triple-negative breast cancer. Breast Cancer Res, 5(12):1-15.
- Chalise P, Fridley BL (2017). Integrative clustering of multi-level omic data based on non-negative matrix factorization algorithm. PLoS One, 12(5):e0176278.
- Tibshirani, R., Walther, G., Hastie, T. (2001). Estimating the number of data clusters via the Gap statistic. J R Stat Soc Series B Stat Methodol, 63(2):411-423.
- Le DT, Uram JN, Wang H, et al. (2015). PD-1 Blockade in tumors with mismatch-repair deficiency. N Engl J Med, 372(26):2509–2520.
- Davoli T, Uno H, Wooten E C, et al. (2017). Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science, 355(6322):261-U75.
- Geeleher P, Cox N J, Huang R S (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol, 15(3):1-12.
- Geeleher P, Cox N, Huang R S (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One, 9(9):e107468.
- Rand W M. (1971). Objective criteria for the evaluation of clustering methods. J Am Stat Assoc, 66(336):846-850.
- Vinh, N. X., Epps, J., Bailey, J. (2009). Information theoretic measures for clusterings comparison. Proceedings of the 26th Annual International Conference on Machine Learning - ICML ’09. p. 1.
- Jaccard Distance (Jaccard Index, Jaccard Similarity Coefficient). Dictionary of Bioinformatics and Computational Biology.
- Fowlkes E B, Mallows C L. (1983). A method for comparing two hierarchical clusterings. J Am Stat Assoc, 78(383):553-569.
- Hoshida, Y. (2010). Nearest Template Prediction: A Single-Sample-Based Flexible Class Prediction with Confidence Assessment. PLoS One, 5(11):e15543.